#读取台站和地震信息
import pandas as pd
import numpy as np
import math

eqinfo      = "/home/pierce/project/predict4ts/速报目录.xls"


def selposition4eqid(eqid):
    lat = np.nan
    lon = np.nan
    
    stationinfo = "/home/pierce/project/predict4ts/Table4Dataprocess.xlsx"
    stationposi = "/home/pierce/project/predict4ts/gravitylist.dat"
    df4stinfo   = pd.read_excel(stationinfo, sheet_name='unitlist4js4gPhone')
    df4stposi   = pd.read_csv(stationposi, header=None, sep='\\s+')
    rs4eqid     = df4stinfo[df4stinfo["仪器编号"]== eqid]
    stid        = rs4eqid["台站代码"]
    #print(stid)
    ptid        = rs4eqid["测点代码"]
    #print(ptid)
    if len(stid) >0 and len(ptid) >0:
       #print(stid.iat[0])
       #print("{}".format(ptid.iat[0])) excel 的东西
       rs4stid  = df4stposi[(df4stposi[2] == stid.iat[0]) & (df4stposi[3]=="{}".format(ptid.iat[0]))]
       #print(rs4stid)
       if len(rs4stid) > 0:
          lat = rs4stid.iloc[0,0]
          lon = rs4stid.iloc[0,1]
          #print(rs4stid.iloc[0,0])
          #print(rs4stid.iloc[0,1])
    return lat, lon

def haversine(lat1, lon1, lat2, lon2):
    """
    计算地球上两点之间的大圆距离（震中距）
    :param lat1: 地震震中的纬度（单位：度）
    :param lon1: 地震震中的经度（单位：度）
    :param lat2: 观测点的纬度（单位：度）
    :param lon2: 观测点的经度（单位：度）
    :return: 两点之间的距离（单位：公里）
    """
    # 将经纬度从度转换为弧度
    lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])

    # 计算差值
    dlat = lat2 - lat1
    dlon = lon2 - lon1

    # Haversine 公式
    a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

    # 地球半径（单位：公里）
    R = 6371.0

    # 计算距离
    distance = R * c
    #print(distance) 
    return distance

def seleqinfo(eqid, oblat, oblon, dist4sel, rank4sel, eqinfo, eq4sel):
    #eqinfo      = "/home/mw/input/nas6133/mnt/data/data/weijin/input/zl_spfx/速报目录.xls"
    #eq4sel      = "/home/mw/input/nas6133/mnt/data/data/weijin/input/zl_spfx/{}_EQ.csv".format(eqid)

    #print(eq4sel)
    df4eqcon    = pd.read_excel(eqinfo, sheet_name='速报目录')
    #oblat       = 41.82
    #oblon       = 86.17
    #dist4sel    = 400
    #rank4sel    = 5.0
    nLabel      = '震中距(km)'
    df4seleqcon = pd.DataFrame(columns=df4eqcon.columns.tolist() + [nLabel])
    #print(df4seleqcon)
    for index, row in df4eqcon.iterrows():
        #print(row)
        eqlat    = row['纬度(°)']
        eqlon    = row['经度(°)']
        eqrank   = row['震级(M)']
        dist4eq  = haversine(eqlat,eqlon,oblat,oblon)
        #print(dist4sel)
        #print(dist4eq)
        #print(dist4eq - dist4sel)
        #print(rank4sel - eqrank)
        #地震时间要在我数据的时间范围
        if dist4eq <= dist4sel and eqrank >= rank4sel:
            row[nLabel] = "{:6.2f}".format(dist4eq)
            #print(row)
            #print(row.to_frame().T)
            df4seleqcon = pd.concat([df4seleqcon, row.to_frame().T], ignore_index = True)
            #break;

    df4seleqcon.to_csv(eq4sel,index=False)
    print('保存成功!')

def seleqinfo4somecol(inFile):
    #print(inFile)
    df4eqcon = pd.read_csv(inFile)
    columns_to_extract = ['序号','发震日期（北京时间）', '震级(M)', '震中距(km)', '震中位置']
    NO      = df4eqcon[columns_to_extract[0]].to_numpy()
    EQDATE  = df4eqcon[columns_to_extract[1]].to_numpy()
    EQRANK  = df4eqcon[columns_to_extract[2]].to_numpy()
    EQDIS   = df4eqcon[columns_to_extract[3]].to_numpy()
    EQPOS   = df4eqcon[columns_to_extract[4]].to_numpy()
    return NO, EQDATE, EQRANK, EQDIS, EQPOS

if __name__ == "__main__":
   eqid = 'X212MGPH0096'
   dist4sel = 400
   rank4sel = 5.5
   eqinfo      = "/home/pierce/project/predict4ts/速报目录.xls"
   eq4sel      = "/home/pierce/project/predict4ts/{}_EQ.csv".format(eqid)
   oblat, oblon = selposition4eqid(eqid)
   seleqinfo(eqid, oblat, oblon, dist4sel, rank4sel, eqinfo, eq4sel)

   






